Global synchronization of bursting neurons in clustered networks 
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We investigate the collective dynamics ol bursting neurons on clustered network. The clustered 
network is composed of subnetworks each presenting a small-world property, and in a given sub- 
network each neuron has a probability to be connected to the other subnetworks. We give bounds 
for the critical coupling strength to obtain global burst synchronization in terms of the network 
structure, i.e., intracluster and intercluster probabilities connections. As the heterogeneity in the 
network is reduced the network global synchronization is improved. We show that the transitions 
to global synchrony may be abrupt or smooth depending on the intercluster probability. 



PACS numbers: 05.45.Xt,87.19.1j 

I. INTRODUCTION 

The human brain is a complex network consisting of 
approximately 10^^ neurons, linked together by 10^** to 
10"'^^ connections, amounting to 10^ synapses per neuron 
[l[ . There are evidences that synchronization may be re- 
lated to a series of processes in the brain. Synchronized 
rhythms have been observed in electroencephalograph 
recordings of electrical brain activity, and are thought 
to be an irnportant mechanism for neural information 
processing [2]. Experimental observations reveals syn- 
chronized oscillations of neuron networks in response to 
sensory stimuli in a variety of brain areas Q . Such 
synchronized rhythms reflect the hierarchical organiza- 
tion of the brain, and occur over a wide range of both, 
spatial and temporal scales 

Complex network such as the brain are typically com- 
posed by subnetworks that are sparsely linked, leading to 
the appearance of clusters. Typically, models for brain 
activity deals with clustered network with small-world 
architecture Q. 

Certain neurons in the brain exhibit burst activity: 
bursts of multiple spikes followed by a rest state hyperpo- 
larization. These bursting neurons are important in dif- 
ferent aspects of brain function such as movement control 
and cognition 0-@- 

If the neurons possess two distinct time-scales such as 
spiking and bursting, the bursts (slow time-scale) tend 
to synchronize at smaller synaptic strengths [l3, Ei • We 
consider two or more neurons to be bursting in a syn- 
chronized manner if they start a given burst at nearly 
the same time, even though the fast spiking may not 
be synchronized. We want to analyze the collective dy- 
namics the a clustered network of bursting neurons and 
give bounds for the smallest coupling strength needed to 
achieve global synchronization on the bursting scale. 

We investigate bursting dynamics on excitatory clus- 
tered networks. Close to the global synchronized state 
the phase dynamics of the burst may be described, to first 



order approximation, by a Kuramoto-like model. Further 
investigations in the Kuramoto model allows us to deter- 
mine the critical coupling strength to obtain burst global 
synchronization in terms of the probabilities of intraclus- 
ter and intercluster connections. 

We show that the route to a global synchronized state 
reveals to distinct types of transitions as a function of the 
interclusters probabilities. If the intercluster probability 
is small enough then the transition to global synchro- 
nization is smooth, as opposed to higher values of the 
intercluster probability. In the latter case the network 
may present an abrupt transition to global synchroniza- 
tion. 



II. CLUSTERED NETWORKS 

A small-world network has typically an average dis- 
tance between nodes comparable to the value it would 
take on a purely random network, while retaining an 
appreciable degree of clustering, as in regular networks. 
Typical brain networks also display a small-world prop- 
erty @. 

We generate small-world networks following the pro- 
cedure proposed by Newman and Watts 12], inserting 
randomly chosen shortcuts in a regular network [13,]. To 
verify that the coupled map network built according to 
procedure proposed by Newman and Watts has the prop- 
erties of a small- world network, we have computed the so- 
called clustering coefficient, defined as the average frac- 
tion of pairs of neighbors of a site that happen also to 
be neighbors of each other, and the average separation 
between sites, that they are showed in Fig. ^ by circles 
and squares, respectively, as a function of the probability. 
The average network distance between sites must be of 
the same order as for a random graph, on the other hand, 
the clustering coefficient must be much greater than for 
a random graph. Then, in virtue of these requirements. 
Fig. ^ suggests that the small-world property is dis- 
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played when the value of the probability is around 10 




FIG. 1. Normalized clustering coefficient (o) and average sep- 
aration between node (□) versus probability of non-local con- 
nections, with N = 100 and p* = 10~^. 

We consider a network composed of M small-world 
subnetworks with L neurons each (Fig. [5]). A given node 
has a probability to connect with other nodes inside 
the cluster, and a probability po to be connected with 
neurons outside the cluster, the intercluster probability. 




FIG. 2. Schematic representation of M = 4 one-dimensional 
small- world networks, each of L = 8 nodes. 

Neuron model: There are a number of mathematical 
models which emulate neuronal activity, ranging from 
differential equations [M] to discrete-time maps [15|] . We 
choose the Rulkov model 

x{n + l) = f{x{n),y{n)) = —^—^+y[n), (1) 

y{n + 1) = g{x{n), y{n)) = y{n) - ax{n) - /3, (2) 

where x{n) is the fast and y{n) is the slow dynamical vari- 
able. Furthermore, a affects directly the spiking time- 
scale and is chosen so that the time series of presents 



an irregular sequence of spikes. The parameters a and /5 
describe the slow time-scale. 

Phase Dynamics: We consider the neurons to be non- 
identical, that is, they possess a mismatch in their inter- 
nal parameters, we choose a as a mismatch parameter. 
Complete synchronization in networks of non-identical 
neurons is not possible. However, weaker synchroniza- 
tion such as phase synchronization [l6| may take place. 

To study phase synchronization, we need to introduce 
a phase for the chaotic dynamics of the bursts. This turns 
out to be a non straightforward task. For chaotic oscilla- 
tors the phase cannot be defined unambiguous. Different 
ways to introduce a phase are possible, each one being 
chosen according to the particular case studied [13, [ill ■ 

We consider the phase being increased by 27r at the 
beginning of each burst. We consider the burst to be- 
gin when the slow variable y{n), which presents nearly 
regular saw-teeth oscillations, has a local maximum, in 
well-defined instants of time we call tk- The duration 
of the chaotic burst, tk+i — i/t, depends on the variable 
x{n) and fluctuates in an irregular fashion as long as x{n) 
undergoes chaotic evolution jljl]. 

We can define a phase describing the time evolution 
within each burst, varying linearly from to t^+i 

(j3{t) ^ 2^k + 2Tr /~*\ ,tk<t< tk+i (3) 

tk+l — Ik 

where tk denotes the time occurrence of the fcth burst. 
The bursts occur in a coherent manner, that is, the time 
interval between burst have a small deviation. 

Chemical Synapses: We analyze a neuron network cou- 
pled via excitatory synapse. We model the synapse as a 
static sigmoidal nonlinear input-output function with a 
threshold and a saturation parameters 20] . We order the 
neurons within the network from 1,2,..., ML, and de- 
note the dynamical variables of the ith neuron as Xi{n) 
and yiin). 

The coupling function is given by 

V{xi{n),Xj{n)) = {xi{n) - Vs)S{xj{n)), 

where (the reversal potential) Vs > Xi{n) for any Xi{n) 
implies that the synapse is excitatory. Here Vg — 2.0. 
The synaptic coupling function S is modeled by the sig- 
moidal function 

^(-) = i + e-M.-e.) - (4) 

We take the saturation parameter A = 10, and the 
value of 9s = —0.25 is chosen in such a way that every 
spike within a single neuron burst can reach the thresh- 
old. 

The network dynamics is described by 

N 

x.,{n+ 1) = f{xi(n),yi(n)) + e^aijV{xi{n),Xj(n)), 
yi{n+l) = g{xi{n),y,{n)), (5) 
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where i = 1,2,..., ML, a = [3 = 0.001, and ai to be 
different for each node on uniform values deviate in the 
interval [4.1,4.4]. e > is the overall coupling strength, 
and A — (a^ ) is the adjacency matrix, aij is 1 if neuron 
i is connected to neuron j, and zero otherwise. 

The local part of the interaction considers the nearest 
and the next-to-the-nearest neighbors of a given node. 
Moreover, with probability pi a new long range connec- 
tion is created inside the each cluster, and with proba- 
bility Po to create an outside connection. 



III. GLOBAL PHASE SYNCHRONIZATION OF 
BURSTING NEURONS 

The collective behavior of the network may be cap- 
tured by the network mean field. 

^ ML 

1=1 

The emergence of a collective behavior enhances the 
mean field dynamics. 

To illustrate the mean field dependence on the synchro- 
nization we choose M = 4 and L = 100 with pi = 10~^, 
Po = 10~^. For small values of the overall coupling pa- 
rameter e = 0.01 the neurons burst in a nonsynchronized 
fashion. The mean field presents no dynamics, but a 
noise-like signal, fluctuating around the mean value Fig. 

The more synchronized neurons the higher the mean 
field. For a coupling e = 0.04 the mean field already 
captures the dynamics of the burst (slow time-scale) dy- 
namics Fig. [Hljb). For the chosen pi and po, if the cou- 
pling strength is large enough then the network globally 
burst phase synchronize. Close to global synchroniza- 
tion e = 0.1 the mean field already reveals the bursting 
dynamics, while the fast oscillation due to the spike are 
filtered, see Fig. [2Jc). This shows that even though the 
burst are synchronized the spikes can remain out of syn- 
chrony. 

The order parameter characterizes the emergence of 
synchrony in the network 
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-^expO^j)), 
j=i 



(6) 



where i?„ and $„ are the amplitude and angle, respec- 
tively, of a centroid phase vector for a one-dimensional 
network with periodic boundary conditions. If the burst- 
ing phases (f>n are spatially uncorrelated, their contribu- 
tion to the result of the summation in Eq. (jH]) is small. 
However, in a globally phase synchronized state the order 
parameter magnitude asymptotes the unity. 

The time averaged order parameter magnitude is given 

by 
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FIG. 3. Behavior of the mean field as a function of the overall 
coupling parameter with fixed pi — 10~^, po ~ 10""^, M — A 
and L = 100. a) e = 0.01 mean field presents no dynam- 
ics, but a noise-like behavior, b) e = 0.04 some neurons are 
synchronized and the mean field exhibit slows oscillators, c) 
£ = 0.1 close to global synchronization the mean field shows 
the bursting dynamics. 



clearly depends on e and both pi, po in and out con- 
nection probability. If the bursting dynamics is globally 
synchronized then Rk, 1. 

Our numerical analysis shows that clustered networks 
with good global synchronization properties is possible 
only for when a Po is larger than a certain critical value. 
This critical values clearly depends on the number of sub- 
networks. For a network that contains many clusters a 
larger po is required to achieve global synchronization. 

Fig. m shows the evolution of the average global order 
parameter R for fixed pi = 10~^ and L = 100. In Fig. 0] 
a) we have M = 2 and b) M = 4, where we depict the 
behavior of R for various values of the Po- In the absence 
of outside connections and e large enough i? « 1/Af, and 
if Po and e are large enough i? « 1. For e = it must 
yield 'R oc (ML)"^/^ so if L > 1 or M > 1 E must be 
close to zero. 

The global synchronization is achieved when R — I. 
However, if phase synchronization is achieved 4>i ~ 4>j 
an almost global synchronized behavior can be obtain by 
setting a threshold value near 1, we choose R — 0.95. 
The minimum coupling strength necessary to achieve 
R — 0.95 is the critical coupling strength Ec, the critical 
coupling parameter for global burst phase synchroniza- 
tion. 

Phase Reduction: the application of the Kuramoto 
phase reduction techniques [2l|, leads to an approxi- 
mate phase equation to describe the collective phase dy- 
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FIG. 4. Time averaged order parameter as a function of the 
coupling strength for pi = 10~^ and L = 100 fixed, a) M = 2 
and b) M = 4. In both figures we have po = (o), p^ = 
1.6 X 10"* (□), po = 2x 10-3 (A). 



the probabilities of intracluster and intercluster connec- 
tions. The critical coupling strength reads 



' p,{L - I) + po{N - LY ^' 

where the constant C depends on the threshold value of 
the average global order parameter R for defining global 
synchronization, that can be regarded as a fitting pa- 
rameter. Besides, a node is connected with L — 1 nodes 
in the same cluster with probability pi that in our case 
is L — 5, and the node is connected with N — L nodes 
from different cluster with probability po- Thus, up to 
first order approximations, the onset of global bursting 
synchronization can be described in terms of Eq. ([7]) 

Our detailed numerical analysis corroborates the the- 
oretical prediction. In Fig. [S]pi = 10"^ and L = 100, 
furthermore we use M = 2 in Fig. a) and M = 4 Fig. 
b) . In both figures the circles represent the critical values 
in accordance with the threshold global order parameter 
equal 0.95. The solid curve is the theoretical prediction 
of the equation ([7]) with Sc + S, 



namics of the neuron ensemble. Since the burst dynam- 
ics is coherent. It is possible to introduce a coordinate 
change such that the phase of the ith reads 

where ipi depends on the particular neuron and C is a 
small noise term p7j . 

Now, if the neurons are close to the onset of global 
synchronization « c/ji, the phase are described by "0^. 
Assuming that phase gradient is approximately constant 
along the trajectories, the phase reduction up to first 
order approximation renders 
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where uJi depends on the neuron parameters and on the 
node degree. Close to the onset of synchronization ipj — 
ipi w sin('0j — 'ipi)- So as a result the phase dynamics of 
the neuron network can be described (up to first order 
approximation) by the Kuramoto model 

0i = w., + e ^ aij sm{4)j - 
i=i 

where, e is a rescaled overall coupling strength and is 
the adjacency matrix. 

Recent results [l^ have given bounds for the critical 
coupling strength required for global synchronization to 



FIG. 5. £c versus Po for a network with L — 100 and pi = 
10"^ The circles are for the threshold value of the global 
order parameter equal 0.95 and solid line Sc ~ C/[pi{L — 5) -I- 
Po{N - L)] + 5, where (a) M = 2, C = 0.07 and 5 = 0.045, 
and (b) M = 4, C = 0.112 and 5 = 0.032. 

The probability Po increases and the critical coupling 
strength required for global synchronization decreases. 
This shows that as the network become are homogeneous 
the global synchronization properties are enhanced. 
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IV. TRANSITION TO PHASE SYNCHRONY 

The transition to synchronization occurs in an inter- 
mittent fashion, where R has laminar phases. To analyze 
the intermittent behavior we introduce the quantity 

1 

where Nji is the number of occurrences of -R > ^ = 
0.95) that we find within a time interval A„ and s is the 
number of configurations of the connections distributed 
with the determined probability po- Thus, F can be in- 
terpreted as the fraction of global synchronization in a 
given time interval. 

F depends on the coupling strength e, and may present 
distinct behaviors as a function of the intercluster prob- 
ability Po- In Fig. [Sl^a) we depict the dependence of F 
on Po- The transition to global synchronization is abrupt 
for Po = 10~^ and Sc ~ 0.05, while it is not abrupt for 
Po = 10-2 and Po = 2 X lO'^. 

Such behavior can be understood heuristically in terms 
of the synchronization properties of each isolated net- 
work. If Po is small enough, each network is almost in- 
dependent. Hence, each network tend to synchronize by 
itself. So, as we increase the coupling each cluster syn- 
chronize and only then a behavior towards global syn- 
chronization is seen. This leads to a smooth transition. 
However, if po is large then the collective dynamics of 
each subnetwork is strongly affected by the remaining 
subnetwork. So the synchronization occurs in one sub- 
network is quickly spreads over the whole network, lead- 
ing to a abrupt transition. 



by the Rulkov map, furthermore the connective architec- 
ture presents the small- world property. We demonstrated 
the possibility of obtaining bursting global synchroniza- 
tion of Rulkov neurons in a clustered networks, where 
the probability of intercluster is varied and the clusters 
exhibit small- world property. Moreover, the network be- 
comes more synchronizable when the probability of in- 
tercluster links is increased. That is, for a fixed pi the 
increasing of Po^ the out connection probability, enhances 
global synchronization. This means that when hetero- 
geneity between the communities (clusters) is suppressed 
the network as a whole synchronizes better. The critical 
coupling parameter required for global synchronization 
as a function of the probability of intercluster presents 
the same behavior as Kuramoto-type dynamics. 
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FIG. 6. Fraction of time where the clustered network has 
larger R than 0.95 versus coupling strength. We consider a 
clustered network with L = 100, M = 2, pi = 10"^, s = 20, 
and Po = IQ-^ (o) and po = 10^^ (□). This two probabilities 
leads to distinct transitions to global synchronization. 



V. CONCLUSIONS 

In conclusion, we presented a neuronal clustered net- 
work model in which the neuron dynamical is described 
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